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Abstract 

An anisotropic ice flow line model is applied to a flow line through Swiss Camp (69.57 N, 49.28 
W) in West Greenland to estimate the dates of internal layers detected by Radio-Echo Sounding 
measurements. The effect of an anisotropic ice fabric on ice flow is incorporated into the steady 
state flow line model. The stress-strain rate relationship for anisotropic ice is characterized by an 
enhancement factor based on the laboratory observations of ice deformation under combined 
compression and shear stresses. By using present-day data of accumulation rate, surface 
temperature, surface elevation and ice thickness along the flow line as model inputs, a very close 
agreement is found between the isochrones generated from the model and the observed internal 
layers with confirmed dates. The results indicate that this part of Greenland ice sheet is primarily in 


steady state. 
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(1). Introduction 

Radio-Echo Sounding (RES) has been used to investigate subsurface properties (internal structure) 
of ice over widespread areas of the Antarctic and Greenland ice sheets. Internal layering from RES 
measurements (Figure 1) provides a vertical profile in an ice sheet, which could previously be 
obtained only from drilled boreholes. The radar echoes arising from internal layers in ice are 
caused by the sudden changes in complex dielectric permittivity of layers in the ice sheets (Fujita 
and others, 1999). The mechanisms causing the sudden changes have been examined in detail, and 
changes in shallow layers have been attributed primarily to ice density variations. Changes in 
deeper layers have mainly been assigned to changes in electrical conductivity due to acidic fallout 
from volcanic eruptions and/or changes in impurity concentration associated with climatic 
transitions (Clough, 1977; Fujita and Mae, 1994; Fujita and others, 1999; Hammer, 1980; Harrison, 
1973; Millar, 1981 ; Moore, 1988; Paren and Robin, 1975). Variation of ice crystal orientation can 
also cause the changes in permittivity at the large depth of the ice sheet (Harrison, 1973; Fujita and 
Mae, 1994; Fujita and others, 1999). 

RES-detected internal layers are widely recognized as isochrones: former ice sheet surfaces that 
have been buried and deformed by ice flow (Clough, 1977; Dahl-Jensen and others, 1997; Hammer, 
1980; Morse and others, 1998; Nereson and others, 2000; Whillans, 1976). Therefore they would 
contain information about the ice sheet history from which inferences can be made about some of 
the climate changes. For example, by knowing the ages of the internal layers the past 
accumulations were inferred from the horizontal variation in thickness between the layers on Taylor 
Dome, East Antarctica (Morse and others, 1998) and around Summit of Greenland (Dahl-Jensen 
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and others, 1997; Fahnestock and others, in press). With the lack of the time information, Nereson 
and others (2000) estimated the pattern of accumulation over Siple Dome, West Antarctica, from 
observed pattern of internal layering. 

To date the internal layers, the most common method is to combine them with the annual-layer 
thickness measurements from the associated ice cores. A system has been recently developed 
(Fahnestock and others, in press) to trace the internal layers from RES imageries and then to extend 
the ages of the layers at the GRIP ice core site at the Summit of the Greenland ice sheet, based upon 
the available GRIP age-depth relationship (Johnsen and others, 1997), to the same layers along the 
airborne flight lines. The dated internal layers mainly cover the Northern areas of the Greenland ice 
sheet (Abdalati, in preparation), which provides a good reference to validate the numerical model 
calculation. In some regions where the internal layers can not be detected or can not be extended 
from the dated ice cores, such as, ablation zone, numerical modeling may give a good 
approximation to the ice age at any depth of ice sheets. 

In this paper, we apply the anisotropic ice flow line model developed by Wang and Warner (1999) 
to a flow line through Swiss Camp (69.57 N, 49.28 W) in West Greenland to estimate the dates of 
internal layers obtained by Radio-Echo Sounding technique (Fahnestock and others, in press). The 
available measurements of dated internal layers are used to validate the model calculation. 

(2). The Model 

The model used here is an anisotropic steady-state ice flow line model (Wang and Warner, 1999). 
The stress-strain rate relationship for anisotropic ice is characterized by an enhancement factor, 
defined as the ratio of tertiary strain rate for anisotropic ice to minimum strain rate for isotropic ice, 
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based on the laboratory observations of ice deformation under combined compression and shear 
stresses (Li and others, 1996). The model has been described in detail by Wang and Warner 
(1999). Here we give a brief review about the flow relations used in the model to describe 
anisotropic ice flow. 

In this two-dimensional model, the horizontal transverse flow is neglected so that the flow is 
entirely constrained to the vertical compression and horizontal shear along the flow line in steady 
state balance. If x denotes the direction of flow and z the vertical direction the flow relations 
between strain rate and stress for components of shear and compression are 

e„=Gi\K(nrl r„ (1) 

e z = F(\)A 0 {T)t 2 0 t : , (2) 

where t xz is shear strain rate, s, compressive strain rate, A 0 (T) temperature dependent parameter, T 
temperature, r 0 octahedral shear stress, t xz shear stress, r. ’ compressive deviatoric stress, G(A. C ) and 
F(A. C ) are enhancement factors of shear and compression component, respectively, and a c is a 
compression factor defined as 

^=r;/(ri+r;T 2 • ( 3 ) 

Laboratory experiments indicate that G(A. C ) and F(?, c ) can be simplified to be equal to E(^ c ) (Li and 
others, 1996) according to: 

E(A C ) = E, (E c / E s ) K , E s = 1 0 and E c =3, (4) 

E s and E c are respective enhancement factors for shear or compression alone. While \c varies from 
1 to 0 as the stress situation varies from purely confined compression (shear stress Xxz=0) to simple 
shear (compressive deviatoric stress r^=0) enhancement factor E(X C ) increases from 3 to 10. 
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In equations (1) and (2), the temperature -dependent coefficient A 0 (T) is based on the laboratory 
experiment results (Budd and Jacka, 1989) tabulated by Wang and Warner (1998, Table 2) which 
has the similar values calculated from Arrhenius relation (Paterson, 1994). The shear stress is taken 
as the driving stress 

Txz=PgaZ, (5) 

in terms of the ice density p, the acceleration due to gravity g, the depth Z and the surface slope a. 
The octahedral shear stress is taken as 


r B =[j <ri + T])t 


( 6 ) 


on the assumption of the ice flow corresponding to a confined vertical compression stress combined 
with a horizontal shear stress. 


After rearranging above equations, a cubic equation for the shear strain rate, s x: , is finally obtained 
as 

4 -f E(A c )A 0 (T)t 3 J£ 2 xz + £:]) = 0 . (7) 

This equation involves vertical compressive strain rate, £ z , which can be calculated from horizontal 
velocity based on the assumption that ice is incompressible. 

Equation 7 is the stress-strain rate relation for anisotropic ice flow used in the flow line model and 
solved iteratively for the shear strain rate (see detail in Wang and Warner, 1999). 


(3). Application of the model to the flow line 

Flow lines over whole Greenland ice sheet were generated based on the assumption of the ice 
flowing down slope in the direction perpendicular to the surface contours using 5 km gridded 
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surface elevations derived from European remote-sensing satellite (ERS-1) radar altimeter data 
(Zwally and Brenner, 2001). The flow line studied in this paper (Figure 2), from the ice divide 
ridge through the Swiss Camp Station towards northwest to the coast, was chosen because an 
airborne radar survey run very close to the flow line. The internal layers were traced from RES 
imagery (Figure 1) along the flight line and dated by extending the same layers from GRIP ice core 
site where dating was available, which were used to validate the model calculation. Furthermore, 
the previous filed observations, such as surface velocity (Hofmann, 1975), accumulation (Benson, 
1962) and temperature (Mock and Weeks, 1966), were available near this flow line and provided 
the input data to the model. 

The Swiss Camp (69.57 N, 49.28 W) located at the West Central Greenland was established near 
the equilibrium line in 1990, estimated from surface balance measurements of the 1980's. The 
surface mass balance has been studied by means of the GPS measurements. The modeling study 
along the flow line through the Camp provides information about the internal dynamics of ice in the 
vicinity of the Swiss Camp. 

Figure 3 shows the data profiles used as inputs to the flow line model. Swiss Camp Station is 510 
km from the ice divide along the direction of flow. 

The ice sheet topography along the studied flow line is shown in Figure 3a. Surface elevations used 
in the model were taken from 5 km grids satellite radar altimeter data (Zwally and Brenner, 2001) 
and bedrock elevation by subtracting ice thickness (Bamber and others, in press) from the surface 
elevation. Surface and bedrock profiles obtained from RES measurements along the flight line are 
plotted in Figure 3a in the section with the observed internal layers. Comparisons with the RES 
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measurements, show the irregular basal topography, but the profile from the 5 km grids appears 
reasonably to represent the smoothed topography incorporated into the model. 

12 internal layers are traced from RES imagery (Figure 1) and displayed in Figure 3a. The 
discontinuous internal layers between 410 and 430 km from the ice divide are due to unclear RES 
imagery in this section, but the pattern is then resumed beyond 430 km. The corresponding ages of 
the layers shown, obtained by tracing the continuous layers back to the dated GRIP core site, vary 
from 2,339 to 12,329 years as presented in Table 1 . 

Internal layers along the flow line seem to reflect the shape of the bedrock undulation with 
decreasing amplitude as the distance from the bedrock increases. This can be seen more clearly in 
Figure 1 . Near the Greenland Summit along the ice-divide ridge, however, where the bedrock is 
smooth, the internal oscillations are caused by variations of the dynamic velocity fields rather than 
the bedrock undulations (Dahl-Jensen and others, 1997). 

Surface accumulation rate along the flow line used in the model is a trend line shown in Figure 3b. 
This was determined based on the previous observations (Benson, 1962) at the locations near the 
flow line (see Figure 2) and 50 km grids database (Zwally and Giovinetto, 2000). The 
accumulation rate is zero at Swiss Camp station on the equilibrium line and is negative at the lower 
elevations where ablation exceeds precipitation. 

The surface temperature profile, shown by solid line in Figure 3c, was interpolated from a 50 km 
grids database (Giovinetto, personal communication, 2001), showing a good agreement with the 
previous observations (Mock and Weeks, 1966) at the locations near the flow line (see Figure 2). 
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The surface velocity profile (Figure 3d) is assumed to be the observations along the Expeditions 
Glaciologiques Internationales au Groenland (EGIG) traverse line (Hofmann, 1975), which runs 
close to the flow line studied in this paper (see Figure 2), based on the assumption that surface 
horizontal velocities along the flow lines vary slowly along the surface elevation contour lines. 
Using laboratory-based flow relations results in an overestimate of surface velocity (~8 times 
higher) by integrating shear strain rate down to the bedrock. This is because the model does not 
consider the reduction of enhancement and shear stress near the bedrock. This reduction may cause 
the reduced shear strain rates near the bedrock which has been found from several of the borehole 
inclination measurements in Antarctica (e.g. Russell-Head and Budd, 1979; Etheridge, 1989; 
Morgan and others, 1998). In Greenland it has been found at Dye 3 borehole (Dahl-Jensen, 1985; 
Dahl-Jensen and Gundestrup, 1 987) that the shear strain rates reach the maximum at the bottom due 
to a high concentration of dust and other impurities, but at the depth between 200 m and 25 m 
above the bed the shear strain rates are almost constant along with the reduction of the enhancement 
factors. To compensate for the overestimated velocity without turning any parameter in the flow 
law, we adopted a reasonably simple scheme from Wang and Warner (1999) by terminating 
integration of shear strain rate near the bedrock to match the measured surface velocity. 


Basal sliding velocity was estimated based on the study of the ice sheet modeling along a flow line 
jointed with the central flow line of Jakobshavn Isbras drainage basin (see Figure 2) (Funk and 
others, 1994) using Equation (Huybrechts and De Wolde, 1999) of 


sliding 



( 8 ) 


where a constant A s =1.8xl0 10 N' 3 yr" 1 m 8 , tb is basal shear stress and Z* the height above 
buoyancy. 
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All input data were interpolated so that the horizontal resolution was I km. A rescaled vertical 
coordinate was used by subdividing the ice thickness into 100 evenly spaced bands. The basal 
temperature gradient of 0.02 °C m' 1 (Funk and others, 1994; Kostecka and Whillans, 1988) was 
used in the calculation of the steady state temperature. 

(4). Results and Discussion 

(a). Comparison of observed isochrones with observed internal layers 

The ages of the ice at any depth of the ice sheet along the flow line were determined by numerically 
integrating the velocity fields following the ice particle trajectories from the surface, using the 
present-day surface accumulation rates and assuming a stead state ice flow. 

Several isochrones calculated from the model are displayed in Figure 4 to compare with the 
observed internal layers. Using the smoothed bedrock and a trend surface accumulation rate results 
in smoothed modeled isochrones, don't capture the high-frequency oscillations of the observed 
internal layers, but agree quite well on larger scales. 

The model simulating isotropic ice flow has been run by replacing the flow relation for anisotropic 
ice with the relation for isotropic ice, i.e. using Equation 7 with E(X. C )= 1 in the model calculations. 
Under all same conditions the isochrones calculated from the model for isotropic ice flow do not 
match the observed internal layers well, which can be seen in Figure 5 and Figure 6a. 

Figure 5 gives a quantitative expression of the ages computed from the models for both anisotropic 
and isotropic ice with the ages obtained from GRIP cores for 12 observed internal layers. The 
modeled ages at the depths of the observed horizons are calculated as a function of .v and shown in 
the bars. The mean values are obtained by averaging those modeled ages along each observed 
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layer. The larger error bars in lower layers are due to the effect of the bedrock undulation, which 
will be substantially reduced by using flight line bedrock (rough curve in Figure 4) in the models. 
Figure 5 clearly indicates the importance of including anisotropy in ice flow modeling. 

(b) . Relation of age and depth, and reconstructed surface accumulation rates 

The comparison of the relations of age and depth derived from the model results with the 
observations is shown in Figure 6a for a point 200 km from the ridge. All three profiles in Figure 
6a closely agree in the approximate top 1/3 ice thickness. At further depth, the modeled profile for 
anisotropic ice still follows the observations closely in contrast with the obvious departure shown 
by the dashed line for isotropic ice, which indicates that isotropic ice is no longer appropriate when 
anisotropy of the ice crystal fabric develops with the increasing of the depth and that it is important 
to consider the effect of the ice crystal anisotropy on the modeling of ice flow in the deeper parts of 
ice sheets. 

Based on the assumption of steady state ice flow, we reconstructed the accumulation rates for the 
internal layer nearest to the surface of the ice sheet with 2,339 years. The age-depth relation given 
in Figure 6a and the modeled vertical strain were used in the calculation. The reconstructed 
accumulation rates are found close to the present-day values as shown in Figure 6b, confirming that 
this part of the Greenland ice sheet has not experienced major changes in accumulation rate and the 
ice flow status. 

(c) . Ice flow at Swiss Camp Station at/near the equilibrium line 
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In this study, present-day surface accumulation rate, surface temperature, surface elevation, ice 
thickness and surface velocity were used as model inputs. The ice flow regime through the depth 
was modeled along the flow line. The agreement between the modeled isochrones and the observed 
internal layers implies that the model incorporated with ice crystal fabric anisotropy has given a 
proper description for the ice flow. Here, we choose Swiss Camp site, at/near the equilibrium line, 
as an example to show the modeling results. The mass balance of the Greenland ice sheet at Swiss 
Camp site will be studied based on these results and the surface observations. 

Vertical profiles of several outputs from the model at Swiss Camp are summarized in Figure 7. The 
parameters shown are: shear and vertical compressive strain rates, horizontal and vertical velocities, 
enhancement, stresses, temperature, and age. 

Shear strain rates ( £ xz ) increase with depth, which is accompanied by an increased enhancement 
(E(^ c )) from 3 to 10, as the range of stress regimes change from predominantly vertical 
compressive stress (rj) in the upper layers to predominantly horizontal shear stress (x xz ) near the 
bed. The transition layer from the compression stress dominance to the shear stress dominance is 
about 300 meters deep (Figure 7g). 

It is important to note several specific points: 1) Since accumulation rate is zero at the site the 
vertical velocities (V z ) represent the total horizontal advections. 2) Temperature profile (T) shows 
that the bottom ice reaches the melting point, indicating a basal sliding is existing. 3) The age-depth 
profile (Figure 7h) shows that the ice near the bottom is exceeding 20,000 years old. 
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(5). Conclusion 

The application of an anisotropic steady-state ice flow line model to the flow line in West 
Greenland is presented in this paper. A close agreement is found between the isochrones generated 
from the model using present-day input data and observed internal layers with confirmed dates, and 
this strong agreement indicates that this part of Greenland ice sheet is essentially in steady state. 
This portion of the ice sheet has been close to its present form for a period at least 12,000 years. 
This success comparison clearly indicates that the model is a good approximation for the flow of 
the ice sheet, and that the ice crystal fabric anisotropy must be considered in order to accurately 
model deeper ice flow. 
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Figure Captions: 

Table 1. Ages of the internal layers displayed in Figure 3a and Figure 4. 

Figure I. RES imagery showing clearly internal layers along the airborne flight line indicated by 
dashed line in Figure 2. 

Figure 2. Map of flow lines in West Greenland showing the contours of surface elevation in 

meters. The studied flow line through the Swiss Camp Station is shown in thick line, the 
airborne flight line in dashed line, and EGIG traverse line provided surface velocity 
observations (Hofmann, 1975) in dotted line. The locations of the measured 
accumulation (Benson, 1962) and temperature (Mock and Weeks, 1966) are indicated by 
solid circle. The sites of GRIP borehole and Swiss Camp Station are marked by asterisk. 

Figure 3. Profiles of input data used in the model showing the location of Swiss Camp Station 

indicated by asterisk, (a). Flow line topography (smoothed curves) used in the model, the 
measured topography from airborne flight (rough curves) and internal layers (thin lines) 
represented the ages in Table 1 . (b). Surface accumulation rates: dots indicate previous 
measurements (Benson, 1962) at the locations shown in Figure 2, the dashed line is 
interpolated from 50 km grids database (Zwally and Giovinetto, 2000), and the solid line 
is the trend of profile used in the model, (c). Surface temperatures: dots indicate 
previous measurements (Mock and Weeks, 1966) at the locations shown in Figure 2, solid 
line is interpolated from 50 km grids database (Giovinetto, personal communication, 

2001) used in the model, (d). Surface horizontal velocities: dots indicate the 
measurements (Hofmann, 1975) along EGIG traverse line shown in Figure 2, asterisk 
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shows the observation at Swiss Camp site and the solid line is the profile used in the 
model. 

Figure 4. Isochrones (solid line) calculated from the model and observed internal layers (dashed 
line) traced from RES imagery 7 (Figure 1). Surface and bedrock elevations used in the 
model are plotted in smoothed curves and the RES observations in rough curves. 

Figure 5. Plot of the ages computed from the models against observed ages from the layer 

measurements. Open triangles and open circles indicate mean values of computed ages 
from the models for anisotropic and isotropic ice along each measured layer, respectively. 
Bars show the range of variations between computed and observed ages along each layer. 
The dashed line indicates the ideal match situation. 

Figure 6. (a). Relation of age and depth, at 200 km from the ice divide, generated from the models 
for anisotropic ice flow (solid line) and using Equation 7 with E(7ic)=l for isotropic ice 
flow (dashed line), and from the observed internal layers (solid circle), (b). Surface 
accumulation rate used in the model (thick line) and reconstructed from the observed 
internal layer nearest the surface of the ice sheet with 2,339 years (thin line). 

Figure 7. Vertical profiles of several outputs from the model at Swiss Camp Station, 510 km from 
the ice divide. They are: (a), shear strain rate ( s x: ), (b). horizontal velocity (V*), (c). 
enhancement factor (E(X C )), (d). temperature (T), (e). vertical compressive strain rate 
( e , ), (0- vertical velocity (V z ), (g). shear stress (t xz ), compression stress deviator (r.) and 
octahedral shear stress (t 0 ), and (h). age. 
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Table I . 
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POPULAR SUMMARY 
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Modeling of Ice Flow and Internal Layers along a Flow Line through Swiss 
Camp in West Greenland 

W.L. Wang, HJ. Zwally, W. Abdalati, and S. Luo 

Snow falling on the Greenland ice sheet forms layers that sink in time as more 
snow continues to be added. While the ice layers sink, they also move toward the 
edge of the ice sheet as the ice flows under the force of its own weight. The depth, 
age, and spacing of these layers provides valuable information for studying the age 
of the ice sheet, the rates of ice flow, the ice-sheet history, and its past climate. 

The depth of the layers have been measured by low-frequency airborne radars that 
sense the internal radar reflections from the layers, as well as the bedrock at the 
bottom of the ice sheet. A numerical model of ice flow, which includes various 
rates of surface snow accumulation and the ice-flow properties at various 
temperatures and depths, is used to determine the dates of the observed layers 
12,000 years backward in time. A novel feature of the model is the use of 
enhanced rates of ice flow based on the laboratory measurements of ice 
deformation under combined forces of compression and shear. Comparison 
between the modeled and measured layer depths for an ice flowline in West-central 
Greenland indicates that part of the ice sheet has been close to the present 
condition of balance between accumulation and ice flow for about the last 12,000 
years. In the future, the model will be used with new data on snow accumulation, 
temperature, and surface elevation changes from satellites to study how the ice 
sheet may be changing with on-going changes in climate. 


